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Abstract 

0^ ■ 

In the last few years the hydrodynamic formulation of quantum mechanics, equiv- 
^ ■ alent to the Bohmian equations of motion, has been used to obtain numerical solu- 

te . tions of the Schrodinger equation. Problems, however, have been experienced near 

wave function nodes (or low probability regions). Here we attempt to compute wave 
' functions and Bohmian trajectories for the interference of one particle or of two 

CD ■ identical particles. It turns out that the large number of nodes (i.e. interference 

minima) makes the hydrodynamic equations impractical, whereas a more straight- 
forward solution of the Schrodinger equation gives very good results. 
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1 Introduction 



Bohmian trajectories were first proposed in an attempt to restore determinism 
in quantum mechanics [1]. In Bohm's view, quantum particles have at every 
instant well-defined positions which, however, can only be known probabilis- 
tically. The particles follow deterministic trajectories governed by equations 
of motion similar to Newton's, except that a specific quantum contribution 
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must be added to the classical potential. This quantum potential explicitly de- 
pends on the particles' total wave function and is responsible for characteristic 
quantum-mechanical effects like barrier penetration. All statistical predictions 
of quantum mechanics can be obtained through averages over trajectories [2,3]. 

One of the first numerical computations of Bohmian trajectories was carried 
out in the context of two-slit interference [4]. It showed vividly how one-particle 
interference effects can be understood in terms of particle dynamics. Trajec- 
tories associated with two-particle interference were also shown explicitly to 
reproduce standard quantum-mechanical results [5]. 

Bohmian trajectories can usually be computed rather straightforwardly if the 
particles' wave function is known analytically. Computation can also be carried 
out from a numerical approximation to the exact wave function. But this, it 
turns out, can be viewed from a different perspective. In the past few years, the 
computation of Bohmian trajectories has led to a powerful way of numerically 
integrating the Schrodinger equation. 

The method is closely connected with the hydrodynamic formulation of the 
Schrodinger equation, which goes back to the early years of quantum mechan- 
ics [6,7]. In this formulation, the evolution of the wave function is associated 
with that of a fluid whose motion can be obtained through the trajectories 
of its elements. Among the quantum-mechanical problems that have been 
adressed in this way are photodissociation [8,9], reactive scattering with the 
Eckart barrier [10,11], the quartic double-well potential [12], and the harmonic 
oscillator with quartic anharmonicity [13]. The method has a number of ad- 
vantages, in particular the use of a relatively small number of grid points and 
its applicability to higher-dimensional problems. It does, however, have dif- 
ficulties in dealing with regions where the wave function vanishes or nearly 
vanishes. 

The purpose of this paper is to investigate the numerical computation of wave 
functions and Bohmian trajectories in the context of particle interference [14]. 
This specifically quantum-mechanical phenomenon illustrates perhaps more 
than any other the properties of quantum superpositions. Moreover, inter- 
ference minima correspond to zeros or near-zeros of the wave function, and 
therefore make severe tests on numerical methods. The hydrodynamic method 
and algorithms for its numerical solution will be reviewed side by side with 
the method based upon separating the Schrodinger equation into its real and 
imaginary parts. Both will be used to make detailed computations of the wave 
functions and Bohmian trajectories associated with the interference of one 
particle or two identical particles. Comparison with results obtained through 
an exact solution of the Schrodinger equation will show that zeros (or nodes) 
of the wave function make the hydrodynamic computation prohibitive in com- 
puter resources, whereas the approach using the real and imaginary parts of 
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the wave function yields accurate results in reasonable time. 



2 Wave functions and trajectories 



The Schrodinger equation for a system of n particles interacting through a 
potential V is given by 

where is the mass of particle i and the Laplacian operator with respect 
to that particle's coordinates. We will focus here on the interference of one 
particle or of two identical particles. There is then only one mass and, in 
dimensionless units, Eq. (1) can be written as 

^^ = -\^'^ + V{f)^p. (2) 

The vector r now stands for the position in configuration space, and the op- 
erator for the Laplacian in that space. 



2.1 Hydrodynamic equations 



The hydrodynamic equations follow from the Schrodinger equation when the 
wave function is written in polar form. They were first used as a basis for the 
numerical solution of Eq. (1) a number of years ago [15], but the method has 
been substantially improved recently. 



To get the hydrodynamic equations we substitute i/j = \/ P cxp{iS) in Eq. (2), 
where S and ^/P are real dimensionless functions and VP is nonnegative. 
Equating separately the real and imaginary parts of (2), we obtain 

^ = -V-{PVS), (3) 
f = -^(V5f-Q-nr), (4) 

where Q, the quantum potential, is given by 

Q = -^^^- (5) 
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The Bohmian trajectories of the particles are defined by writing the following 
equation for the velocity in configuration space: 

v^VS. (6) 



We now substitute (6) in (3) and (4). Noting that VaVb — V^Va (a and b are 
coordinate indices), we find 

^^-P{t)V-v{t), (7) 
^ = -V(Q + y), (8) 
where the Lagrangian derivative is defined as 

Eqs. (7) and (8), together with 

Df 

must be solved to get the Bohmian trajectories and, eventually, the wave 
function. 



The numerical solution of Eqs. (7), (8), and (10) can be carried out in two 
different ways. The first one, called Lagrange's viewpoint, uses the fact that 
the Lagrangian derivative represents the total time derivative with respect to 
a moving coordinate system. Grid points are then defined which move with 
the particles according to Eq. (10). This method has a number of advantages. 
First, the Bohmian trajectories are automatically calculated. Secondly, as grid 
points follow Bohmian trajectories, they remain concentrated in regions of high 
probability. As a consequence, fewer points are needed since the grid adapts 
itself throughout time evolution. Finally, use of the Lagrangian derivative gives 
differential equations with very few terms. 

Euler's viewpoint, in contrast with Lagrange's, uses a fixed grid. Eqs. (7) 
and (8) are solved with (9) substituted into them. That method may be more 
fiexible for the purpose of computing spatial derivatives. 



2.2 Schrddinger equation 



Several methods for the numerical solution of the Schrodinger equation in 
its original form were developed over the years [16,17,18,19]. Heller's ap- 
proach [18], in particular, was used for the analysis of atom diffraction by 
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surfaces through computation of Bohmian trajectories [20,21]. Here we shall 
separate the Schrodinger equation into its real and imaginary parts [16] by 
substituting i/j = i1^r + iipi in (2). We get 

(n, 

^ = ivVii - Vi,^. (12) 

Starting with the value of ^0 at a given time, these equations are to be solved 
over an interval of time. Bohmian trajectories are then computed using (6) or, 
explicitly, 

'^^"^'^ = = V (arctan f . (13) 



dt ' ' I 



3 Numerical methods 



In this section, we address the problems of appropriately evaluating spatial 
derivatives and carrying out time integrations, in each of the two schemes 
considered. 



3.1 Hydrodynamic equations 



In the Lagrangian viewpoint, one needs to compute spatial derivatives on a 
grid that changes at each time step. We use the moving weighted least squares 
method proposed in Ref. [10]. It consists in fitting a series of polynomials to 
values of functions in a neighborhood of the point were the derivative is to be 
evaluated. Provided that the neighborhood is small enough, the function to be 
differentiated should be well represented by low-order polynomials. Function 
derivatives will then be given by the coefficients of the polynomials. This 
method can be adapted to almost any point distribution. 



To be more specific, suppose we want to evaluate derivatives of a function / 
at a point vq. We first write / as a finite series of polynomials around tq, that 
is. 



M 



To)- 



(14) 



s=l 

We now use the values of / at A^;, neighboring points f„ of fo, and find the 
coefficients Ug by minimizing the following expression: 



n=l 



f{rn) - Tl^iasPsir^ 



n 2 



(15) 
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Introducing the rectangular matrix 

A„. = (16) 

and the vector 

bn = (17) 

we find that is minimized if 

A^-6 = A^-A-a. (18) 

This is called the normal equation [22]. The weighted least squares approxi- 
mation is implemented by solving (18) for the unknown a. Note that • A 
is a square matrix. Once the are known, derivatives of / can be evaluated 
from Eq. (14). The standard error (T„ is determined by assigning larger weights 
to closer points, using for instance a Gaussian distribution around tq. Further 
details on the use of the weighted least square method in connection with the 
hydrodynamic equations can be found in Refs. [11,12]. 

Time integration in Lagrange's viewpoint is based on the following discretiza- 
tion of the Lagrangian derivative: 

Dm fit+At)-fit) 

"dT ^ At ■ ^^^^ 

Here / is evaluated on the moving grid, whose points follow particle trajecto- 
ries. In Euler's viewpoint we have 

dm fit + At) - fit) 

dt At ' ^^"^ 

where / is now evaluated on a fixed grid. 

The spatial derivatives turn out to be smoother if wc make the transformation 
P = exp{2g). Making use of Eqs. (7), (8), and (10), we find in Lagrange's 
viewpoint 

rn{t + At) = fn{t) + AtVnit), (21) 
Vn{t + At)^Vnit)-AtW{Qn + Vn), (22) 

gn{t + At) = gn{t) - ^At V • Vn, (23) 



where 



Qn = -^^^ = -\ [ (V^n)' + ^'9n] ■ (24) 



2 VPn 

Here fn{t), for instance, stands for the value of the function / at the grid point 

— * — » 

n at time t, and Vgn = (Vg)n- 
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In Eulcr's viewpoint, the following equations have to be solved, together 
with (24): 

Vnit + At) = vn{t) - AtViQn + K) - At{v^it) ■ V}Vn{t), (25) 
Qnit + At) = gn{t) -^AtV-Vn- Atvn ■ V^„. (26) 

It should be pointed out that a higher-order scheme like the conventional 
Runge-Kutta method cannot be used for time integration with the weighted 
least squares method, since the functional form of the spatial derivatives 
changes at each time step. 



3.2 Schrddinger equation 



The numerical solution of Eqs. (11) and (12) requires discrete approximations 
to the second-order spatial derivatives of ipR and ipj. Let A denote the grid 
spacing. We use the following approximation, which neglects terms of order 
A^ and higher: 

-{/(x + 2A)+/(a;-2A)}]. (27) 

This formula cannot be used near the grid boundaries, where /(x ± A) and 
f{x ± 2A) may not be defined. In this case we write, for instance, 

= mix) - 15Af{x + A) + 214/(a; + 2A) 

- 156/(x + 3A) + 61f{x + 4A) - 10f{x + 5A)] , (28) 

+ Uf{x + 2A) - 6f{x + 3A) + fix + 4A)] , (29) 
with similar expressions on the other side of the grid. 

A fourth-order expression for first derivatives is also used for the computation 
of Bohmian trajectories. 

Once the spatial discretization is done, Eqs. (11) and (12) read as 

d 1 

Q^'^Rn = -l^Fn{lpIm), (30) 

d 1 

^V'/n = -T^Fnii^Rm)- (31) 
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For a grid with points, this makes up a system of 2N coupled first-order 
differential equations. From Eqs. (27)-(29), one can see that for a given value 
of n, the index m assumes up to six different values. 

Since oscillations of the real and imaginary parts of the wave function may be 

important, the numerical integration of V'i? and ipi requires an accurate and 
stable scheme. We use the fourth-order Runge-Kutta method. Note that we 
improve on Ref. [16] in both the spatial and the time discretization. 



4 One- and two-particle interference 

An idealized interference setup is shown in Fig. 1, where parameters later to 
be used in wave functions are indicated. The source S either emits one particle 
at a time, which may go through one of the slits and be detected on the screen. 
Or it emits two identical correlated particles, with identical x momenta and 
opposite y momenta, so that if one particle goes through slit A the other goes 
through slit B. 

Let and il^siriji) be the partial wave functions for particle i going 

through slit A or B. Just like the symmetry of the setup, we assume that ipA 
and ipB transform into each other under reflection through the x axis, that is, 

i/jAixi, Vi, t) = iJB{xi, -yi, t). (32) 
The z coordinate is omitted throughout. 

For one-particle interference, the global wave function is given by 

*one(r, t)=N' [Mn, t) + Mn, t)] , (33) 

where the configuration space coordinate r corresponds to the one-particle co- 
ordinate fi. Here A/" is a normalization constant. For two-particle interference, 
r— (fi, r2) and we can write 

*two(r, t)^^f [^A(ri, t)V'B(r2, t) ± Mn,t)Mr2, t)] . (34) 

The + sign corresponds to bosons, for which the global wave function is sym- 
metric under particle exchange, while the — sign corresponds to fermions, for 
which the wave function is antisymmetric. In the one-particle case the inter- 
ference pattern shows up on the screen, whereas in the two-particle case it is 
a property of configuration space. 

At t = 0, the partial wave functions are picked as plane waves in the x- 
direction, and Gaussian wave packets in the y-direction, centered on the ap- 
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propriate slit. Explicitly, 



i^A{ri,t^ 0) = (27r(To) ^ exp 



iy^ - Yf 



(35) 



with ipB given through (32). The time evolution of such wave functions in free 
space (where y = 0) is known exactly. It is given by 



^.(n.() = (2.a?)-"%xp|-(?^ + . 



where 



kit 



(36) 



(37) 



For plane waves along x, one can show [5] that the x-coordinate Bohmian tra- 
jectory is simply given by x{t) — x{0) + k^t. In the numerical implementation, 
we therefore concentrate only on the ?/-coordinates. One-particle interference 
thus reduces to a one-dimensional problem, whereas two-particle interference 
is a two-dimensional problem. 

We recall that Bohmian trajectories have been obtained, from exact wave 
functions, for one-particle interference in [4] and for two-particle interference 
in [5]. In the numerical computations of wave functions and trajectories that 
follow, we let throughout F = 1, ctq = 0.2, and k^ = 0.1. 



5 Results and discussion 



5.1 Hydrodynamic equations 



One of the main advantages of the Lagrangian viewpoint is the possibility of 
concentrating grid points in regions of high probability. Accordingly, our first 
attempts at solving the hydrodynamic equations for one-particle interference 
used grid points in the immediate neighborhood of the slits only. The numeri- 
cal results obtained with such initial conditions were very different from what 
should be expected. Instead of building up an interference pattern, they repre- 
sented essentially independent Gaussian wave packets emerging from each slit. 
Clearly then, grid points are needed in the whole region between the slits, even 
where the probability of finding a particle is very low. This is related to the 
fact that the development of plateaux and troughs in the quantum potential 
responsible for the formation of fringes really begins around y = [4]. These 
remarks point to one of the main problem encountered with the hydrodynamic 
equations: wave function nodes [12,13,23,24]. 
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The reason why nodes or quasinodes of the wave function arc apt to cause 
problems in the hydrodynamic approach is apparent from Eq. (5). While the 
denominator of the quantum potential then nearly vanishes, the numerator 
normally does not. The quantum potential may thus experience rapid and 
important variations which challenge approximation procedures. Moreover, in 
the Lagrangian approach the trajectories tend to group in regions of higher 
probability, thereby going away from nodes. Thus precision is lacking just 
where the most important variations of the quantum potential occur. The 
problem is especially acute in our case of one-particle interference. Nodes then 
correspond to some of the most interesting parts of the wave function, namely 
interference minima. 

Since the wave function (33) for one-particle interference nearly vanishes ini- 
tially at the point |/ = 0, it should help to understand the node problem. With 
grid points in the neighborhood of slits only, no interference pattern is formed 
and trajectories go through the node just as if it were absent. But Bohmian 
trajectories should never cross, hence in this case the quantum potential is 
not properly calculated. When points are added near the node, the quantum 
potential can be calculated better. Fig. 2 shows that getting an accurate ap- 
proximation is no easy matter. As expected, the approximation tends to get 
better as points are added and more polynomials are used. In the Lagrangian 
approach, however, the matrix • A has to be inverted at every point and 
every time step. Since the dimension of the matrix is equal to the number 
of polynomials, computation time increases quickly, which makes the method 
inefficient. 

Moreover, small inaccuracies in the initial quantum potential cause important 
effects in later times, as can be seen in Fig. 3. Although a good approximation 
is found for the quantum potential with 801 grid points between —4 and 4 
and fifth-order polynomials, oscillations appear in the quantum potential and 
in the velocity as early as t = 0.01 (1000 time steps). These oscillations can 
be caused by instabilities in the time integration or by small errors in the 
approximation of the spatial derivatives, or both. 

The grid used for these tests was uniform. Although the weighted least square 
method allows for nonuniform grids, it would not really help concentrating 
points around the initial node. Other nodes would develop as the interference 
pattern builds up, which would also require additional points. 

Several solutions have been proposed in connection with the node problem, 
for example grid adaptation [25] and a hybrid method consisting of solving 

the Schrodinger equation near nodes and the hydrodynamic equations else- 
where [24,26]. In the case of interference, however, nodes are permanent in 
time as well as moving in space, so that adaptation is too expensive. As far as 
the Schrodinger equation goes, we shall show in the next section that it is in 
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fact more efficient than the hydrodynamic equations, and therefore does not 
need to be coupled with it. 

Instead of using Lagrange's approach, we may try Euler's. Since grid points 
are then fixed, precision can be controlled easily. Because matrix inversion 

is carried out only at the ffist time step, computation times are somewhat 
shorter. Yet use of the weighted least square method still makes this approach 
inefficient. As a matter of fact, the behavior of the quantum potential and 
velocity is roughly similar in the Eulerian approach as in Figs. 2 and 3, since 
almost no dispersion of the wave packet has yet occurred at t = 0.01. In 
addition, numerical instabilities arise due to the terms in Eqs. (25) and (26) 
that are absent in the Lagrangian approach. This suggests that a higher-order 
scheme is probably needed for time propagation. 

Other ways of approximating derivatives can be used with Euler's approach 
and give good results in reasonable time. The approximations described in 
Sect. 3.2, for example, give very good results for the initial node. However, 
they require more points than needed in connection with the Schrodinger 
equation, and the time integration is highly unstable. In this case the fourth- 
order Runge-Kutta scheme could be used. But again the Schrodinger equation 
looks more promising, since it involves at most second-order (instead of third- 
order) derivatives. 

There is no point here to look at two-particle interference, since the behavior of 
V, Q, and P is similar to what was found for one particle, and the two-particle 
problem requires a much larger number of points. 

5.2 Schrodinger equation 

Figs. 4 and 5 show the real and imaginary parts of the wave function for 
one-particle interference. Excellent agreement is found between the numerical 
solution of the Schrodinger equation (dots) and the exact value (solid lines). 
The grid spans y- values between —13 and 13 with spacing A equal to 0.1, for 
a total of 261 points. 5000 steps were used to go from t — to t = 1. 

With the real and imaginary parts of the wave function in hand, Bohmian 
trajectories can be obtained through Eq. (13). Some of these trajectories are 
shown in dotted lines in Fig. 6, solid lines representing trajectories obtained 
from the exact wave functions. Again the agreement is excellent. In general 
the grid should cover enough space to avoid boundary effects. But with the 
scheme of Sect. 3.2, second-order derivatives are computed quickly and, even 
with 5000 times steps, the overall computation time (for the wave function and 
trajectories) is a few minutes on a 1.8 GHz Pentium 4 processor. For compar- 
ison, the computation time of wave functions in the Lagrangian approach is 
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between one and two orders of magnitude higher [27] . 

As mentioned earher, two-particle interference here is a two-dimensional prob- 
lem, and therefore the grid must be much larger. We have used a square grid 
with 261 X 261 points, both ^-coordinates going from —13 to 13. The com- 
putation time is therefore much longer. Bosons wave functions (the -|- sign in 
Eq. (34)) were used throughout. 

Once again results obtained for the real and imaginary parts of the wave func- 
tion are in very good agreement with exact values. Some Bohmian trajectories, 
obtained with 15 000 time steps, are shown in Fig. 7 as dotted hues. Solid lines 
represent trajectories computed with exact wave functions. The agreement is 
usually very good, with just a small difference showing up in the upper curve 
in Fig. 7a. The main reason for this is here again that the trajectory goes 
through a node of the wave function, where the velocity varies quickly. The 
complete numerical solution for the wave function and trajectories is then 
more demanding. Similar behavior was observed for one-particle interference, 
but in both cases differences are small, in sharp contrast with results from the 
hydrodynamic equations. 



6 Conclusion 



To our knowledge, Bohmian trajectories in the context of interference through 
wave packet spreading have hitherto been calculated only from exact wave 
functions. In this paper, we have investigated two different methods for the 
numerical computation of wave functions. The hydrodynamic equations, writ- 
ten cither in Lagrange's or Euler's viewpoint, arc sensitive to the evaluation of 
derivatives of the quantum potential. This is especially delicate near wave func- 
tion nodes, inevitable in interference problems. Lagrange's viewpoint is usually 
attractive because the grid automatically adapts to regions of high probability, 
and because it addresses higher-dimensional problems with relative ease. Yet 
here proper evaluation of derivatives with the weighted least squares method 
requires a large number of points and high-order polynomials, and time propa- 
gation tends to be unstable. In Euler's viewpoint, derivatives can be computed 
more efficiently, and higher-order schemes can be used for time propagation. 
But then it is simpler and more accurate to use the Schrodinger equation di- 
rectly, which involves at most second-order derivatives (instead of the gradient 
of the quantum potential). In both one- and two-particle interference prob- 
lems, the direct numerical solution of the Schrodinger equation thus provides 
an accurate and relatively quick way of obtaining wave functions and Bohmian 
trajectories. 
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Figure captions 



Figure 1. Two-slit interferometer. 

Figure 2. Initial quantum potential, (a) 12 neighbors, 401 points, order of 
polynomials varied; (b) 12 neighbors, fifth-order polynomials, number of points 
varied. 

Figure 3. Velocity at t = 0.01. Other conditions same as in Fig. 2. 
Figure 4. Real part of the wave function at t = 1. 
Figure 5. Imaginary part of the wave function at t — 1. 

Figure 6. Bohmian trajectories in one-particle interference. The x-coordinate 
is proportional to the ^-coordinate, since x{t) — k^t. 

Figure 7. Bohmian trajectories associated with pairs of particles, (a) |/i(0) = 1, 
^2(0) = -0.6; (b) yi(0) = 1, 7/2(0) = -1.4. 
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Figure 2 
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Figure 3 
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Figure 6 
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Figure 7 
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